# -*- coding: utf-8 -*-
"""
Created on Thu Mar 23 09:28:13 2023

@author: lenovo
"""

from osgeo import gdal
import os
# import glob
import shutil
from tqdm import tqdm # 运行计时


town=['商洛市','商州区','洛南县','丹凤县','商南县','山阳县','镇安县','柞水县',]


vartype=['水源涵养','土壤保持','洪水调蓄','空气净化','水质净化','固碳价值','释氧价值',
         '气候调节','负氧离子','储量碳价值','水资源','landtype']



shapedir = r'F:\公司项目\陕西\商洛\shp'

tiffdir = r'F:\公司项目\陕西\商洛\data'


outpath = tiffdir + '\\masked'

if os.path.exists(outpath):
    shutil.rmtree(outpath)
    os.mkdir(outpath)
else:
    os.mkdir(outpath)
    

for ci in tqdm(town):
    for tp in vartype:
        input_raster=tiffdir+'\\merged\\'+tp+'.tif'
        input_shape=shapedir+'\\'+ci+'.shp'
        output_raster=outpath+'\\'+tp+'_'+ci+'.tif'
        
        
        ds = gdal.Warp(output_raster,
                      input_raster,
                      format = 'GTiff',
                      cutlineDSName = input_shape,      
                      cropToCutline = True,
                      copyMetadata=True,
                      creationOptions=['COMPRESS=LZW',"TILED=True"],
                      dstNodata = -9999)
        ds = None

            





